VANDERMONDE FACTORIZATIONS OF A REGULAR HANKEL 
MATRIX AND THEIR APPLICATION ON THE COMPUTATION OF 

BEZIER CURVES 

LICIO H. BEZERRA* 

Abstract. In this paper, a new method to compute a Bezier curve of degree n = 2m — 1 is 
introduced, here formulated as a Bcrnstcin-Hankcl form in C m , that is, each coordinate of the curve 
is of the form e^-BJ^ (s)HB'f n (s) T e m , where BJ^(s) is a m X m lower triangular Bernstein matrix 
and H is a Hankel matrix. The method depends on Vandcrmonde factorizations of a regular Hankcl 
matrix, and so we begin with a proof, which utilizes Pascal matrices techniques, that given a regular 
Hankel matrix H, there is a finite set of complex numbers 7 such that x m — Pm-l£ m_1 — ••• ~P0 has 
multiple roots, where (po ...p m —i) = (h m +i ... h n H —1 . Therefore, a Vandermonde factorization 
of H can be accomplished by taking a complex number at random, and the Bernstein-Hankel form 
can be easily calculated, thus yielding points on the Bezier curve. We also see that even when H is 
nearly singular, the method still works by shifting the skew-diagonal of H. By comparing this new 
method with a Pascal matrix method and Casteljau's, we see that the results suggest that this new 
method is very effective with regard to accuracy and time of computation for various values of n. 
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1. Introduction. Let H be a Hankel matrix of order n, i.e., (Vi, j G {1, n}) 
Hij = hi + j-\. A very known theorem says that, if H is nonsingular, then a Van- 
dermonde matrix V and a diagonal matrix D exist such that H = VDV T . There 
is a proof of this fact in [5], which utilizes a class of matrices arisen in the theory 
of root separation of algebraic polynomials, namely the class of Bezoutians. Here, in 
section $ 2, from a procedure that is currently utilized in linear prediction to estimate 
parameters in exponential modeling, it is showed that the spectrum of the companion 
matrix C — C{x 1 ), where a; 7 is the solution of the linear prediction system Hx — y~ fl 
with y 1 — (h n +i ••• h^n-x l) T , is simple for all but a finite set of 7. For the values 
belonging to this finite set, there is a more general factorization: H — V C DV^ , where 
V c is a confluent Vandermonde matrix and D is a block diagonal matrix, as it can be 
seen in [3] . Our approach to the proof of the Vandermonde factorization of a nonsin- 
gular Hankel matrix is very similar to the one found in [7], but the proofs are distinct. 
For instance, we make here use of generalized Pascal matrices to quickly obtain some 
general properties of Hankel matrices. 

In section $3, we see that a Bezier curve of degree n — 1, where n = 2m — 1, 
can be described as a Bernstein-Hankel form on C m . Also, in this section a new 
algorithm to compute Bezier curves is proposed, from a Vandermonde factorization 
of the associated Hankel matrix. In section $ 4, results of numerical experiments are 
presented, which strongly suggest that we can compute those curves in a very fast and 
precise way. That is corroborated from the comparisons done with the Casteljau's 
method ([5]) with various values of n. On the other hand, however, several experi- 
ments indicate that the computation of Vandermonde factorization of a Hankel matrix 
is sensitive to its condition with respect to inversion. However, once its skew-diagonal 
entries are shifted toward skew-diagonal dominance the precision of the computation 
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improves, which is a simple and efficient way to deal with the instability of Vander- 
monde factorization of ill-conditioned Hankel matrices, at least for the computation 
of Bezier curves from this approach. 

2. Vandermonde factorizations of a nonsingular Hankel matrix. 

/ hi h 2 ... h n ^ 



Let H 



\ 



Suppose H is nonsingular. Let x 7 be the 



1n—l h n ... h 2n -2 
h n h n+ i ... h 2n -l 

solution of the linear prediction system Hx = y 7 , where y 1 = (h n+ i ... h 2n -i l) T ■ We 
want to show that the set of 7 G C for which the companion matrix C 7 = compar^x^) 
is not diagonalizable is finite. S is a nonderogatory matrix, it suffices to show 

that S, the set of scalars 7 such that the spectrum of C 7 is not simple, is finite. This 
means that, out of this set, the characteristic polynomial of C 7 , p 7 (x), doesn't have 
multiple roots. If a = (ao ... a„_i) T and b = (bo ... 6 rl -i) T are the respective solutions 



of Hx = e„ = (0...01) J and Hx = (h n+1 ... h 2n -i 0) J 



r(x) — r ys(x), 

where r(x) — x n — & n _ix ,l_1 — ... — b\x — bo and s(x) — a„_ia; ,l_1 + ... + ayx + ao- It 
is not difficult to see that S is finite iff r(x) and s(x) don't have any common root. 



then p 1 (x) 
-1 1 



LEMMA 2.1. Let H be anxn nonsingular Hankel matrix. If a = (ao...a„_i) T and 
b = (bo-..b n -i) T are the respective solutions of Hx — e n and Hx — (h n+ i ... h 2n ^i 0) T , 
then ao 7^ or bo 7^ 0. 

Proof. 

Suppose \H(1 : n — 1,2 : n)| 7^ 0. Therefore, from Cramer's rule, ao 7^ 0. Let 
Xi, ...,x n -i be the unique scalars such that 



Xl 



( h 2 



V h„ 



%7l—l 



I K \ 




( K + i \ 


\ ^2n-2 J 


- 


\ h 2n -i ) 



Hence, x — (xo xi ...x„_i) T = 7a + b is the solution of Hx = (h n+ i ... h 2n ~i 7) T , 
with xo = 0, iff 7 = xih n+ \ + ... + x n -ih 2n -i. For other complex numbers 7, 
xo = 7^0 + bo ^ 0, that is, ao 7^ or bo 7^ 0. Notice that ao 7^ 0, and b = iff 
xih n+ i + ... + x„-ih 2n -i = 0. 

Now, suppose H(l : n — 1, 2 : n) = H(2 : n, 1 : n — 1) is singular. First, since H 
is nonsingular, the dimension of span{H(2 : n, 1), H(2 : n,n — 1),H(2 : n,n)} is 
(n — 1), as well as the dimension of span{H(l : n — 1, 1), H (1 : n — 1, n — 1), H(l : 
n— 1, n)}. Hence, iJ(2 : n, n) ^ span{H(2 : n, 1), iJ(2 : n, n — 1)}, whose dimension 
is n — 2. On the other side, H(2 : n, n) e span{H(l : n — 1, 1), i?(l : n — 1, n)} = 
span{H(l : n — 1, 1), i/(2 : n, 1), i/(2 : n,n — 1)}, and so, there exist xo, x„„i, 
where Xo is different from zero and unique, such that 

/ "n+i \ / hi \ ( h2 \ ( hn 

: = x ; + Tl ^ + ••• + ; 

\ ^2n-l / \ / \ h n J \ h 2n _ 2 

Observe that, in this case, for all 7 £ C, xq = 6q 7^ 0, and ao = 0. □ 



From the above proof, there can be at most one complex number 7 such that 
p 7 (0) = —bo — jao = We can also conclude from the lemma [2TT1 that zero is not a 
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common root of r(x) = x n — b n _ix r ' 



/ hi 



Now, define H" — 



b\X — b and s(x) = a„_ 1 .x" 1 + ... + 



. Since H is nonsingular, 



h n ■■■ h 2n -i 7 
\ h n+ i ... 7 k J 
is also nonsingular iff n ^ kq = (h n+ i ... h 2n -i 7) H~ 1 (h n +i ... h 2n ^\ l) T , which is 
equal to (h n+1 ... h 2n -i j)(bo + 7 a o ••• &n-2+7 a ™-2 K-i + ia n -i) T ■ 



Note that H* 



\ 



-bo - ja ^ 








= (k-Kq) 





bn-l - 7 fl n-l 




1 J 




V 1 / 



Therefore, lemma [2TT1 can 



be rewritten as the following lemma: 



Lemma 2.2. Let H" 



h n h n +i \ 



K 
\ h n+ i 



«2n-l 

7 



7 



fee a Hankel matrix, 



where H = H*(l : n, 1 : n) is nonsingular. Suppose that is also nonsingular, that 
is, k ^ (h n+ i ... ft,2n-i 7) H~ l {h n +i ■■■ hin-i l) T ■ Let p be the solution of H*x = 
e„+i. Then, except for one possible complex number 7, po ^ 0. 

Now, let a be any complex number and q~f(x) = p 7 (x + a) = r(x + a) — 7s(x + a). 
In an analogous way to the proof for a — 0, it will be shown that r(a) and s(a) cannot 
be both null because there can be only one complex number 7 such that g 7 (0) = 0. 
To prove this, we introduce some notations and definitions in the following. 

Definition 2.3. Let a G C. P„[ck] fee the n x n is the lower triangular matrix 
defined for each i,j G {1, 2, . . . , n} by 



a 



, otherwise. 



P„[a] is said a generalized lower triangular Pascal matrix. If a = 1, P n [l] = P n is 
called the n x n lower triangular Pascal matrix. 

Some results about these matrices (see [BJ, pQ) are listed in the following lemma: 

Lemma 2.4. Let P n [a] a generalized lower triangular Pascal matrix. Then, 

(a) P„[0] =/„; 

(b) P n [a]P n [p] = P n [a + p]; 

(c) {P n [a])- X = P n [-a]; 

(d) Let a ^ and let G n {a) be the n x n diagonal matrix such that, for all 



k G {!,..., n}, (G n (a)) kk = 



„fc-i 



Then P n [a] = G n (a)P n G n (a) 1 = 
G„(a)P„G„(a _1 ). In particular, P" 1 = G n (-l)P n G n (-l). 

Definition 2.5. For s G [0, 1], the n x n Bernstein matrix B^(s) is the matrix 
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defined for each i,j 6 {1,2, ... ,n} as follows: 

lB e (s)] .._( (jliy-Hl-*)^ , fori>j; 
^ S >^-\ , otherwise. 

A very important fact about Bernstein matrices, which will be used here later, is the 
following proposition, whose proof can be found in pQ: 

Proposition 2.6. Let s € [0, 1] and let B e (s) be anxn Bernstein matrix Then, 
B^(s) = P n G n {s)P~ x , where P n is the n x n lower triangular Pascal matrix and 
G n (s) = diag([l,s,..., S n - 1 ]). 

In the following, we present some relations between Pascal and Hankel matrices. 

Lemma 2.7. Let H be a n x n Hankel matrix and let P n be the n x n lower 
triangular Pascal matrix. Then P n HP^ is still a Hankel matrix. 

Proof. The lemma obviously holds when n = 1. Suppose it holds for all Hankel 
matrices H of order n > 1. Now, let H be a (n+1) x (n+1) Hankel matrix and consider 
P n+1 HPT +1 . Since P n+1 HP^ +1 is symmetric and P n+1 HP^ +1 = [P n HP% v;v T k], 
for some d e C", by induction it suffices to show that, for all A; S {1, ...,n— 1}, 
(P n+ iHpT) n+hk = (P n+ i#Pj +1 ) n , fe+ i. Now, 



J 



rp rp fk - V 

(Pn+iH P n+1 ) n +i,k = e n+1 P n +i ( . )He j+ i = 

3=0 

fc-1 / \ /, ,\ n+k-l 



i=0 j=0 \ / \ J 7 s= 2 1=0 v 7 v 7 

which is equal, from Vandermonde convolution ([5])) to 

s=2 i=0 v 7 v 7 to j=0 v 7 VJ/ 



fe 

i=o 







3+1 



(P n+ iHpT) n , 



,fc+i- 



Corollary 2.8. Let H be a n x n Hankel matrix and a be a complex number. 
Then, P n [a]H P n [a] T is still a Hankel matrix. 

Proof. For a = 0, the result follows from lemma [2~7l Let a ^ 0. Since from lemma 
I2.4I P,, [a] — G n (a)P n G n (a~~ 1 ), where G n (a) = diag (1, a, a n ~ 1 ), it suffices to show 
that G(a)HG(a) is a Hankel matrix. But this is obviously true, for (G(a)HG(a)) i j = 
hi+j-xa^- 2 . □ 

Next we give a proof that r(x) and s(x) don't have any common root by using a 
generalized Pascal matrix technique. 

Proposition 2.9. Let H be a n x n nonsingular Hankel matrix. Let a = 
(ao di ... a„_i) T and b — (bo bi ... 6„_i) T be the solutions of Ha — e n and Hb = 
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(h n+ i ... h 2n -i 0) T , respectively. Then r(x) = x n — b n _ix n 1 — ... — b\X — b and 
s(x) = a„_!a; n_1 + ... + ciix + a don't have any common root. 

Proof. Let 7 £ C and let p 1 = (—b — ja ... — — 7a n _i 1) . Let q 1 — 
(q ... q n -\ 1) T be the vector of coefficients of the polynomial r(x + a) — r ys(x + a). 
We note that g 7 = F n+ i[o;] T p 7 = P n +i [a] T (i? 7 : ) _1 e n+ i, for k = 1 + k . Thus, 
H^P n+ i[a]^ T q 7 = e„+i, and so, 

Hfa = P„ + iH _1 ^^ l+ iH _T 97 = P n +i[-u\H*P n+1 [-a] T q^ = e n+1 . 

is also nonsingular and, from corollary 12.81 is a Hankel matrix. Since iJ 7 = 
Ho+j(e n+ ie^ + e n e^ +1 ) + Ke n+ ie^ +1 , we see that H* = -ff°+7 (e«+ie^ + e„e^ +1 ) + 
(k — 2 n a) e„+ie^ +1 . That is, 

/ ^1 ■•■ h n+ i \ 

/l„ ... /l2n-l 7 
V -7 « / 



H« = 



where -ff 7 (l : n, 1 : n) = H = P„ [— a]-ff-Pn [— a] T is nonsingular and 7 = 7+(iJg)„ + i jr 
Thus, from lemma [2~2l except for one possible complex number 7, (g 7 )o 7^ 0. □ 

Note that (<z 7 )o = only when s(a) 7^ 0, that is, when i? 7 (l : n — 1,2 : n) 
H(l : n — 1, 2 : n) 7^ 0. In this case, 7 = r(a)/s(a). 

Proposition 2.10. Let 7 e C. Let p 7 (a;) = x n - b n _ x x n ~ l - ... - b - 
7(a n _ia;™ _1 + ...+ao) = r(x)— r ys(x) the characteristic polynomial ofC 1 = H\(^)H , 
where -Hi (7) is the Hankel matrix defined by Hi(^/)ek — Hek+i for k = 1, ...,n — 1 
and Hi(j)e n = (h n +i ■■■ h 2n ~i l) T , that is, C 7 = [e^; ej; (/i„+i ... /i 2n -i 7)-/? _1 ] . 
Then the set of scalars 7 such that C 7 is not diagonalizable is finite. 

Proof. C-y is a companion matrix, and hence, a nonderogatory matrix. Thus, it 
suffices to show that the set of scalars 7 such that the spectrum of C 7 is not simple 
is finite. 

Let a E C be an eigenvalue of C 7 , that is, a root of Pj(x). Therefore, r(a) = 
r ys(a). Then, from proposition 12.91 s(a) 7^ 0. So, there are two cases: 

(i) r(a) = 0, and this occurs iff 7 = 0. In this case, Co is not diagonalizable iff 
r'(a) = 0. 

(ii) r(a) 7^ 0, which means that 7 = r(a) / 's(a). Therefore, p 7 (ct) = iff 
r'(cf) = s'(a) — 0, or s'(a) 7^ and r'(ct) = js'(a). 

Therefore, since s ^ and r/s is not a constant, a is contained in the set of the 
roots of r' s — rs', which has at most 2(n — 1) elements. Hence, we can conclude that 
{7 E C I C-y is not diagonalizable} is finite and has at most 2(n — 1) elements. □ 

We can now state the following theorem: 

Theorem 2.11. Let H be a n x n nonsingular Hankel matrix. Let r(x) = 
x n — b n -ix n ~ 1 — ... — bo and s(x) = a n _ia;" _1 + ... + oq, where a = (ao a\ ... a n -i) T 
and b — (bo b\ ... 6„_i) T are such that Ha = e n and Hb = (h n+ \ ... hm-i 0) T - Let 
S = {a E C\(rs' - r's)(a) = and s(a) ^ 0} and T = {r(a)/s{a)\a E S}. 
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Then, for all 7 S C — T, H = V^D^VT , where V 1 = vander{a.\, ct n )i -C7 = 
diag(V~ Hex), { a%, a n } = A(C 7 ), and C 7 is the companion matrix whose last 
row is (b + ja ... 6„_i + 7<2 n -i)- 

Proof. From proposition 12.101 for all 7 6 C — T, A(C 7 ) is simple. Suppose 
{ ai, a„ } = A(C 7 ). Let w = (h n+ i ... hm-i l) T and Hi = [H(2 : n, :); w]. Then, 

C 7 = H 1 H~ 1 = V 1 diag([a u a n ]) V' 1 , 

where V^a — (l a, ... a" _1 ) T , for all i £ {1, n}. So, 

V- 1 ^ = diag([a 1 ,...,a n ])V- 1 H. 
Let d = (di ... d n ) T — V~ x He\. Hence, for all i £ {1, n}, 

V~ x Hei = diag{[a x , ...,a n ]f~ l d = (dxoff 1 ...d n a£, X ) T = D^V^a. 

□ 

3. Bezier curve as a Hankel form. Efficient methods to compute Bezier 
curves of degree n — 1 ([3]) are fundamentals tools in Computed- Aided Geometric 
Design area. The Casteljau's algorithm is a widespread method for this computation. 
However, for each s 6 (0, 1) it demands 0(n 2 ) multiplications. For n not very large, 
there are more efficient methods, like the ones introduced in |10j . where the compu- 
tation of points on these curves is carried out by generalized Ball curves, or the ones 
presented in [SJ, which use fast Pascal matrix- multiplication. Here we show that we 
can describe a Bezier curve as a Hankel form and, hence, we see that we can eas- 
ily compute points of the curve from a Vandermonde factorization of the associated 
Hankel matrix. 

Let Q = (x ,y ), Qi = (xi,yi), Qn-i = (x n -i,y n -i) be n points in M 2 . 
Bezier has his name on the curve B defined from these n points as follows: 

^)=(^:])=E( n 7 1 ) si(i - srl " Qi ' se[o ' i] - 

Let x — (xq ... £„_i) T and x = (yo ... y n -i) T ■ Then, for each s £ [0, 1], 

h(s) = elB e n {s)x and b 2 (s) = e^B^(s)y, 

where B^{s) is a n x n Bernstein matrix. Thus, from lemma [2~31 for each s £ [0, 1], 

his) = elP n G n {s)P- 1 x and b 2 {s) = elP n G n {s)P- l y. (3.1) 

In the following, we discuss different approaches that make use of (|3.ip to compute a 
Bezier curve. 

We can notice that, if B(s) = Bq Q 1 __,Q n _ 1 (s) denotes the Bezier curve deter- 
mined by the points Qq, Qi, Q n -i, then 

B(s) = {l-s)B Qo Q 1 ...Q n _ a (s) + sB Ql Q a ...Q n _ :L (s). 
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Without loss of generality, from now on we will suppose that n, the number of control 
points of a Bezier curve, is odd: n — 2m — 1, m > 1 . In this case, it is easy to conclude 
by induction that, for all k = 0, m — 1, 



k <k 



3=0 

Particularly, for k = m — 1 we have 

m—l 



m— i / _ I \ 

B w = E ( ,• )( l -*r" 1 "^' B ^^ + i-^ +m - 1 w 

j=0 \ J s 



and so, 



m— 1 / , 

'm—l 



j=0 



m-1 _ i\ 

j=o V J / 

where denote the column vectors (ajj . . . Xj+ m -i) T and 

(j/j . . . j/j+ m -i) , respectively, for j = 0, m — l. However, 

E ( ,• K 1 - *) m - 1 - J Ve^P m G(t)P- 1 a;j -... J - +m _ 1 = 
j=o \ •? / 



&i M = E( j ](i-»r 1 - 3 Vej;p m G(i)p; i Ij ... j+m . 1 , 



elP m G(t)P, 




and X^j=o ("j ) (1 — s ) m 1 ^ s ^ x 3—j+m-i is a column vector whose ith coordinate is 
e^ n P m G(t)P~ 1 Xi~i... m+ i-2- Thus, we can state the following lemma, from which we 
can conclude that each coordinate of a Bezier curve is a Hankel form: 

Lemma 3.1. Let n = 2m — 1, where m is an integer greater than 1 and let B{s) = 
(b\(s) &2(s)) T be a Bezier curve of degree n — 1 defined from n points Q = (x~,yo), 
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(x 1 ,y 1 ), Qn-i = {x n -i,y n -i) in Mr. Then 

bi(s) = elB e m (s)H x (B< m ( s )) T e m and b 2 (s) = elB e m (s)H y (B e m (s)f e m , 

where H x — hankel(C x , R x ) and H y = hankel(C v , R y ) are m x m Hankel matrices 
whose first columns are C x = (xo...x m _i) T and C y = (yo...y m -i) T respectively, and 
whose last rows are R x = (x m _i, ...,x n -i) and R y = (y m -i, ...,y n -\ respectively. 

Corollary 3.2. Let n — 2m— 1, where m is an integer greater than 1. Let B be a 
Bezier curve of degree n—1 definedfromn control points, and let x = (xo...x n -i) T and 
V = (yo---y n -i) T be their respective vector of coordinates. Let H x — hankel(C x , R x ) 
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and H y = hankel(C y ,R y ), where C x = (x ...x m _i) T , R x = (x m -i,...,x n -i), C y = 
(yQ---y m -i) T and R y = (y m -i, y n —i)- If H x and H y are nonsingular, then there 
exist complex numbers d± 7 ...,d n , tx,...,t n , d\,...,d n and t\, ...,t n such that 

m tci 

h(s) = J2M^-s + s.t i ) n - 1 andb 2 {s)=J2 c ii( l - s + s -ti) n ~ 1 - ( 3 - 2 ) 



Proof. If H x is nonsingular, from theorem 12.111 there is a Vandermonde matrix 
V = vanderQti, ...,t n ]) and a diagonal matrix D — diag([d±, ...,d n ]) such that H x = 
VDV T . So, 

h{s) = elB e m (s)H x (B e m ( S )fe m = e T m B* m {s)V DV 7 '(B^s)) 7 e m = 



= J^diil - s + s.ti) 2m - 2 = di(l - s + s.tif- 1 , 

i=l i=l 

for e T n B% / {s)Ve i =Y%=a{l-s) m - 1 -i. S i4 = (1 - s + s^)™" 1 for all z 6 {l,...,m}. 
In an analogous way, we conclude that 



b 2 (s) = j2d i (i-s+s.i i ) n - 1 , 



for some d\, d n and t±, t n . D 

The following proposition is about another representation of a Bezier curve of 
degree n — 1 . 



Proposition 3.3. Let n = 2m — 1, where m is an integer greater than 1 and 
let B(s) = (h(s) &2(s)) T be a Bezier curve of degree n — 1 defined from n points 
Qo = (x ,yo), Qi = (xi,Vi), ■ Qn-i = (x n -i,y n -i) ofM 2 . Then 



b i( s ) = X/ afe \ , j sk an d h{s) = ^ h ( " , 1 where 

k=0 ^ ' k=Q ^ 



Oq ■■■ CLn-l 



( b 




- P~ 1 H P- T 

1 m ±1 V ± m ■ 



= P m 1 H x P m 1 and 

\ a m -\ ... a„-i J \ b m -\ 

Proof. Let A = P m 1 H x P m T and B = P- 1 H y p- T . From lemma GHU it follows 

that 

h(s) = e 7 P m G(s)AG(s)P^e m and b 2 {s) = e 7 P m G(s)BG(s)P 7 e m . 
Now, e T m P m G{s) = (p- 1 ) C:!)^ 1 )' Therefore ' 

M,-g.(i(V)(r;)U 
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□ 



We have just proved a property of the Pascal matrix-vector multiplication, which 
is remarked in the following corollary: 

Corollary 3.4. Let n — 2m — 1, where m > 1, let x — (x n ...x n _i) T be a vector 
of C ra and let H x the Hankel matrix defined by {H x )^ — Xi + j-2- Then, a — P n x, 
where a = (ao...a„_i) T is such that (P m H x P^) = a.i+j-2. 

Proof. Let y = G n (-l)x. Let a = P n x = G„(-l)P,7 1 G n (-l)x. Then, 

n-l 

elP n G n {s)P- l y = e T n P n G n {s)G n {-l)a = ^(-l) fe a fe 

k=0 

On the other hand, from proposition [331 (— l)' i+:, ~ 2 ai + j_2 = (Pm 1 HyP r ^ T ) . There- 
fore, cn+j-2 = (P m H x P^). r □ 

4. Numerical experiments. We are going to compare Bezier curves of degree 
(n — 1) computed from the classical Casteljau's algorithm as well as from two other 
descriptions of the curve: as a Hankel form and by using the spectral decomposition 
Bn( s ) = PnG n (s)P~ 1 . We first observe that an uniform scaling of the control points of 
a Bezier curve yields an uniform scaling of the curve and if those points are translated 
by a vector v — (p, q), then the Bezier curve is also translated by v. Hence, without 
loss of generality, we are going to assume that the coordinates of the control points 
are all real positive and also less than or equal to 1. So, we are going to use the 
MATLAB function rand to generate n test control points: A = rand(n, 2). 

The Casteljau's algorithm is a very accurate algorithm to evaluate Bezier curves, 
for it is based on a numerically stable Bernstein matrix-vector multiplication: 



Algorithm 1 Casteljau's algorithm 
n = length(x); 
x = x(:); 
ss = 1 — s; 
for k = 2 : n do 

for t — n : —1 : k do 
x(t) = ss*x(t-l) + s*x(t); 

end for 
end for 
b(s) — x(n) 




10 



L. H. Bczerra 



This multiplication can be seen as a sequence of bi-diagonal matrix- vector multi- 
plications, which becomes well explicit from the following lemma [2]: 
Lemma 4.1. Let B^(t) be a n x n Bernstein matrix. Then 

B°(s) = ££_ x (s)...£^(s) where, for 1 < Jfe < n - 1, 

E%{s) = eiej + ... + e k el + e fe+1 [(l - s)e k + se fe+ i] T + ... + e n [(l - s)e n _i + se n ] T . 

Another way of calculating a Bezier curve is from its description as a Hankel 
form, which allows us to utilize a Vandermonde factorization of the associated Hankel 
matrix, and its algorithm is as follows: We are supposing here that H x and H y 



Algorithm 2 Bezier curve as a Hankel form 

• given n — 2m — 1 distinct points of K 2 , let H x and H v be two mxm Hankel 
matrices formed from their coordinates; 

• choose a number 7 at random and define the vectors x 1 — ■■■ ^7) T 
and y 1 = {h y m+l .../i«7) T ; 

• solve the systems H x z y — x-y and H y w-y — y 7 and consider the companion 
matrices C Zi and C Wi ; 

• find the spectra of C z ^ and C Wj ; 

• define d x = (y x )~ 1 H x e 1 and d y = (Vy)- 1 H v ei, where V x and V» are Van- 
dermonde matrices formed from the spectrum of C Zi and from the spectrum 
C Wy , respectively; 

• for each s £ [0,1], the Bezier curve B(s) is then defined from the equation 



are both nonsingular and that 7 is not one of those numbers which yield in non- 
diagonalizable companion matrices. 

The third way of computing a Bezier curve will be carried out by a Pascal matrix 
method, which computes a Bezier curve B(s) of degree n — 1 via the decomposition 
Bl{s) = P n G n (-s)P n G n (-l): 



Algorithm 3 Pascal matrix algorithm 

• given 71, take t > 1 such that P n (t) is similar to a lower triangular Toeplitz 
matrix T = T(t) with maximum minT^/ maxTij; 

• multiply z = P n G n (— l)x — P n x- and w — P n G n (—l)y = P n y~ via fast 
Toeplitz matrix-vector multiplication; 

• from a Horner-like scheme, evaluate the polynomials e^P n G n (—s)z and 
en P nG n (-s)w. 



We have used a fast Pascal matrix-vector multiplication done from the similar 
Toeplitz matrix T(t) (see [11]). where t has been found by a procedure described 
in [2], plus the B(s) evaluation given by a Horner-like scheme that evaluates the 
polynomial concomitantly with the binomial coefficients. Since the £?(s)-evaluation 
becomes unstable when s approaches to 1, we have introduced a simple procedure to 
improve the evaluation, that is to divide the process of evaluation in two independent 
steps: 

(a) evaluate e^P n G n (—s)z and e^P n G n (—s)w for < s < 1/2; 
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(b) evaluate e^P n G n (—s)z r and e^P n G n (—s)w r for 1/2 > s > 0, which is equiv- 
alent to evaluate e^P n G n {—s)z and e^P n G n (—s)w for 1/2 < s < 1. 

Table 4.1 

Mean run time of computation of 129 points of a Bezier curve of degree N-l by three different 
methods: Casteljau's (C), Hankel form (H) and direct Pascal matrix method (P). The results of the 
second and third methods are compared to the ones obtained by Casteljau's via norm of the difference 
of the computed points by the respective method and by Casteljau's. 



N 


Time (Castcljau) 


Time (Hankel) 


Time (Pascal) 


\\B c -Bh\\ 


\\B C -Bp\\ 


15 


0.005s 


0.007s 


0.004s 


1.3399c- 13 


2.6782c-12 


23 


0.009s 


0.009s 


0.004s 


1.0540c- 11 


2.2427e-09 


31 


0.015s 


0.011s 


0.005s 


2.3082e-09 


1.4962e-06 


39 


0.022s 


0.016s 


0.006s 


9. 7593c- 11 


5.6283e-04 


47 


0.030s 


0.019s 


0.006s 


6.6642e-05 


0.1035 


55 


0.040s 


0.023s 


0.007s 


4.9873e-08 


457.2366 


63 


0.053s 


0.027s 


0.007s 


1.8852e-05 


2.2703e+04 


71 


0.066s 


0.029s 


0.008s 


6.0574e-07 


1.7485c+07 


79 


0.082s 


0.036s 


0.009s 


1.0117c-06 


5.6499e+09 



4.1. Conditioning a Hankel matrix. It is not rare n — 2m — 1 numbers taken 
in the interval [0,1] at random result in an ill-conditioned m X m Hankel matrix H. 
A simple way of handling this is to shift its skew diagonal in order to turn it into a 
skew-diagonal dominant matrix, H = H + <tC, where C is the reciprocal matrix. Let 
Bh and Bjj be the Bezier curves corresponding to H and H , respectively. Then, for 
each s € [0,1], we compute Bh(s) by subtracting a times e^ n B^ n (s)C m (B!^ n (s)) T e m 
from Bfj. Moreover, this quadratic form has a simple formulation as can be seen in 
the next lemma. 

Lemma 4.2. Let C m — hankel(e m ,ej), which is called the reciprocal matrix. 
Then, ifw = e 27ri /"\ 

m 

elB e m (s)C m (B e m (s)) T e m = -^^(l - s + s.w*- 1 )*- 1 . 

Proof. It is easy to see that C m — VDV T , where V = vander(l, w, w™" 1 ) and 
D = diag(l/m,w/m, w 111-1 /m). From the proof of Corollary |3.21 

e T m B e m ( S )VDV T {B e m {s)) T e m = - £ ^(l - s + s^'" 1 )"" 1 - 

□ 

In table 14.21 we can see that this simple technique of preconditioning have im- 
proved the computation of Bezier curves when their control points yield ill-conditioned 
Hankel matrices (cond(H) is the maximum condition number of the two Hankel ma- 
trices formed by the coordinates of the control points). For each Hankel matrix H, 
a was taken as the sum of the absolute values of its entries. Since our Vandermonde 
factorization of a Hankel matrix depends on a value chosen at random, the error 
between the curve computed by Casteljau's and the one computed from that factor- 
ization varied enormously when the Hankel matrices associated with the coordinates 
were ill-conditioned. In table 14.21 for each n, we can see the maximum error among 
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several experiments done. However, sometimes it happened to have a big error fol- 
lowed by a tiny one. Notice that all our experiments have been run in a 32-bits AMD 
Athlon XP 1700+ (1467 MHz). 

Table 4.2 

Mean run time of computation of 129 points of a Bezier curve of degree N-l by three different 
methods: Casteljau's (C), Hankel form (H) and preconditioning Hankel form (PH). The results of 
the second and third methods are compared to the ones obtained by Casteljau's via norm of the 
difference of the computed points by the respective method and by Casteljau's. 



N 


cond(H) 


Time (C) 


Time (H) 


Time (PH) 


||S C -B H ||2 




31 


1.5379c+03 


0.015s 


0.012s 


0.019s 


3.3983c+ll 


2.9510c-ll 


39 


760.4605 


0.022s 


0.016s 


0.024s 


2.6541e+09 


1.1134c-10 


47 


2.9956e+03 


0.030s 


0.019s 


0.031s 


1.3731c+13 


1.0189c-10 


55 


577.1450 


0.041s 


0.023s 


0.036s 


4.3750c+03 


1.7107c-08 


63 


4.2415c+03 


0.053s 


0.027s 


0.042s 


9.8000c+70 


2.5894e-08 


71 


907.6247 


0.066s 


0.029s 


0.049s 


3.7813c+06 


3.2318c-07 


79 


1.1167c+03 


0.081s 


0.033s 


0.057s 


4.1314e+04 


2.1604e-05 



REFERENCES 

[1] L. Aceto AND D. Trigiante, The matrices of Pascal and other greats, Amer. Math. Monthly, 
108 (2001), pp. 232-245. 

[2] L. H. Bezerra and L. K. Sacht, On computing Bezier curves by Pascal matrix methods, 

larXiv:1006l4327l rl [math.NA]. 
[3] P. Bezier, Numerical Control: Mathematics and Applications, John Wiley &: Sons, London, 

1972. 

[4] D. L. Boley, F. T. Luk AND D. Vandervoorde, A General Vandermonde Factorization 
of a Hankel Matrix, in ILAS Symp. on Fast Algorithms for Control, Signals and Image 
Processing, 1997, Winnipeg. 

[5] W. Boehm AND A. Muller, On de Casteljau's algorithm, Comput. Aided Geom. D., 16 (1999), 
pp. 587-605. 

[6] G. S. Call and D. J. Velleman, Pascal's Matrices, Amer. Math. Monthly, 100 (1993), pp. 
372-376. 

[7] M. Fiedler, Special Matrices and Their Applications in Numerical Mathematics, Second ed., 

Dover, Mineola, NY, 2008. 
[8] R. L. Grahan, D. E. Knuth and O. Patashnik, Concrete Mathematics - a Foundation for 

Computer Science, Second ed., Addison- Wesley, Reading, MA, 1994. 
[9] G. Heinig and K. Rost, Algebraic methods for Toeplitz - like matrices and operators, 

Birkhauser, Basel, 1984. 

[10] H. N. Phien and N. Dejdumrong, Efficient algorithms for Bezier curves, Comput. Aided 

Geom. D., 17 (2000), pp. 247-250. 
[11] X. Wang and J. Zhou, A fast eigenvalue algorithm for Pascal Matrices, Appl. Math. Comput., 

183 (2006), pp. 713-716. 



